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Abstract 

Background: Cellular metabolism is one of the most investigated system of biological interactions. While the 
topological nature of individual reactions and pathways in the network is quite well understood there is still a lack 
of comprehension regarding the global functional behavior of the system. To get some insight in this direction, in 
the last few years powerful theoretical methods such as extreme pathways calculation and flux-balance analysis 
(FBA) have been introduced. These methods strongly rely on the hypothesis that the organism maximizes an 
objective function. However only under very specific biological conditions (e.g. maximization of biomass for E-Coli 
in reach nutrient medium) the cell seems to obey such optimization law. A more refined analysis not assuming 
extremization remains an elusive task for large metabolic systems due to algorithmic constraints. 
Results: In this work we propose a novel algorithmic strategy that allows for an efficient characterization of 
the whole set of stable fluxes compatible with the metabolic constraints. Using a technique derived from the 
fields of statistical physics and information theory we designed a stochastic algorithm to estimate the size of the 
affine space containing all possible steady-state flux distributions of metabolic networks. The algorithm, based 
on the well-known Bethe approximation, allows the computation in polynomial time of the volume of a non 
full-dimensional convex polytope in high dimensions. We show results for toy models that compare very well with 
the results of exact algorithms. The result of our algorithm match closely the prediction of Monte Carlo based 
estimations of the flux distributions of the Red Blood Cell metabolic network. It is then used to analyze the 
statistical properties of the average fluxes of the reactions in the E-Coli metabolic network and finally to test the 
effect of gene knock-outs on the size of the solution space of the E-Coli central metabolism. 
Conclusions: We propose a novel efficient algorithmic strategy to estimate the size and shape of the affine space 
of a non full-dimensional convex polytope in high dimensions. The method is shown to obtain, quantitatively and 
qualitatively compatible results compared with known algorithms (where this comparison is possible) being still 
efficient on the analysis of large biological systems, where all other known strategies experience an explosion in 
algorithmic time. 



1 



Background 

Cellular metabolism is a complex biological problem. 
It can be viewed as a chemical engine that trans- 
forms available raw materials into energy or into the 
building blocks needed for the biological function of 
the cells. In more specific terms a metabolic net- 
work is indeed a processing system transforming in- 
put metabolites (nutrients) , into output metabolites 
(amino acids, lipids, sugars etc.) according to very 
strict molecular proportions, often referred as stoi- 
chiometric coefficients of the reactions. 

Although the general topological properties of 
these networks are well characterized, see for exam- 
ple [1—3], and non-trivial pathways are well known 
for many species [4] the cooperative role of these 
pathways is hard to comprehend. In fact, the large 
sizes of these networks, usually containing hundreds 
of metabolites and even more reactions, makes the 
comprehension of the principles that govern their 
global function a challenging task. Therefore, a nec- 
essary step to achieve this goal is the use of mathe- 
matical models and the development of novel statis- 
tical techniques to characterize and simulate these 
networks. 

For example, it is well known that under evo- 
lutionary pressure, prokaryotes cells like E-Coli be- 
have optimizing their growth performance [5]. Flux 
Balance Analysis (FBA) provides a powerful tool to 
predict, from the whole space of phenotipic states, 
which one will these cells acquire. In few words one 
may say that FBA maximizes a linear function (usu- 
ally the growth rate of the cell) subject to biochemi- 
cal and thermodynamic constraints [6] . On the other 
hand, cells with genetically engineered knockouts or 
bacterial strains that were not exposed to evolution 
pressures, need not to optimize their growth. In fact, 
the method of minimization metabolic adjustment 
(MOMA) [7] has shown that knockout metabolic 
fluxes undergo a minimal redistribution with respect 
to the flux configuration of the wild type. Yet, 
in more general situations, the results are unpre- 
dictable, therefore, a tool to characterize the shape 
and volume of the whole space of possible phenotipic 
solutions must be welcome. 

Unfortunately, this characterization has re- 
mained an elusive task. As far as we know the at- 
tempts to obtain such a characterization were always 
based on the Monte Carlo sampling of the steady- 
state flux space [8]. But unfortunately, it appeared 
that this kind of sampling is a very expensive calcu- 
lation and is unsuitable for very large networks. To 



address this problem we propose, using a technique 
derived from the fields of statistical physics and in- 
formation theory, an algorithm that may efficiently 
characterize the whole set of stable fluxes compati- 
bles with the stoichiometric constraints. 



Mathematical Model 

As already mentioned, a metabolic network is an 
engine that converts metabolites into other metabo- 
lites through a series of infra-cellular intermediate 
steps. The fundamental equation characterizing all 
functional states of a reconstructed biochemical re- 
action network is a mass conservation law that im- 
poses simple linear constraints between the incoming 
and outcoming fluxes at any chemical reaction. It is 
call the dynamic mass balance equation: 

f = i + S.,-o (1) 

where p is the vector of the M metabolite concen- 
trations in the network, i (o) is the input (output) 
vector of fluxes, and v are the reaction fluxes gov- 
erned by the M x N stoichiometric linear operator S 
(usually named stoichiometric matrix) encoding the 
coefficient of the M intra-cellular relations among 
the N fluxes. 

As long as just steady-state cellular properties 
are concerned one can assume that a variation in 
the concentration of metabolites in a cell can be ig- 
nored and considered as constant. Therefore in case 
of fixed external conditions one can assume metabo- 
lites (quasi) stationarity and consequently the Ihs of 
[T] can be set to zero. Under these hypotheses the 
problem of finding the metabolic fluxes compatible 
with flux-balance is mathematically described by the 
linear system of equations 

S ■ v = o i = b (2) 

where b is the net metabolite uptake by the cell. 
Without loss of generality we can assume that the 
stoichiometric matrix S has full rows rank, i.e. that 
rank(5) = M, since linearly dependent equations 
can be easily identified and removed. Knowing that 
the number of metabolites M is lower than the 
number of fluxes N the subspace of solutions is a 
(N — M )-dimcnsional manifold embedded in the N- 
dimcnsional space of fluxes. In addition, the pos- 
itivity of fluxes, together with the experimentally 
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accessible values for the maximal fluxes, limit fur- 
ther the space of feasible solutions. This fact may 
be expressed by the following inequalities: 

< v < i/ max (3) 

in such a way that together, [2] and [3j define the 
convex set of all the allowed time-independent phe- 
notipic states of a given metabolic network. 

Sub-dimensional volumes 

Mathematically speaking, the space of feasible solu- 
tions consistent with the equations [2] constitutes an 
affinc space V C R N of dimension N-M. The set of 
inequalities [3] then defines a convex polytope II C V 
that, from the metabolic point of view, may be con- 
sidered as the allowed configuration space for the cell 
states. With the scope of describing it, we will be 
interested in computing the volume of this space and 
certain volumes of subspaces of it. Although concep- 
tually simple, the notion of sub-dimensional volume 
like that of II requires some new definitions. 

Consider any linear parameterization <fi : 
R N - M -> V C R N . A popular choice for (f> is, for 
instance, the inverse of the so called lexicographical 
projection i.e., the projection over the first N — M 
coordinates such that its restriction to V has an in- 
verse. Being </> linear, the (N — M) x N Jacobian 
matrix A is constant and coincides with the matrix of 
4> in the canonical bases. Denoting A = det(ATA)2, 
the Euclidean metric in 1* induces a measure on V 
(which does not depend on <fi): 

J^f{u)du = \J f(c/>(u))du (4) 
allowing to compute the volume of our polytope 
vol 1 /(II) = J l n (v)dv = xj l^-i {n) (u)du (5) 

where In (•) is the indicator function of the set II. It 
is worth pointing out that given the linear structure 
of the metabolic equations, the determinant of the 
mapping is a (scalar) constant. 

Probabilistic framework 

The problem of describing the polytope II can for- 
mulated in a probabilistic framework. We define the 
probability density V as: 

V{y) = voI v (n)- 1 l n ( I /) (6) 



Marginal flux probabilities over a given set of fluxes 
are obtained by integrating out all remaining de- 
grees of freedom. In particular we can define sin- 
gle flux marginal probability densities as integrals 
on the affine subspace W = V fl {v, = ~P} : 



Pi{vi) = 




i.e. the (sub dimensional) volume of the intersec- 
tion between the polytope II and the hyperplane 
{Vi = V}. 

Results and Discussion 
Performance on low dimensional systems 

In this section we will analyze the performance of 
our algorithm against an exact algorithm on low di- 
mensional poly topes. Among the different packages 
available in the Internet, we have chosen LRS [9], a 
program based on the reverse search algorithm pre- 
sented by Avis and Fukuda in [10] that can compute 
the volume of non-full dimensional polytopes. Actu- 
ally, it computes the volume of the lexicographically 
smallest representation of the polytope, that for the 
benchmark used below, coincides with the conven- 
tional volume estimated by our algorithm. 

We have devised a specific benchmark generating 
random diluted stoichiometric matrices at a given 
ratio a — M/N and fixed number of terms different 
from zero K in each of the reactions. All fluxes were 
constrained inside the hypercube < i/$ < 1. As a 
general strategy we have calculated several random 
instances of the problem and measured the volume 
(entropy) of the polytope using the LRS and BP al- 
gorithm. In particular, we have first generated 1000 
realizations of random stoichiometric matrices with 
N = 12, M = 4. Note that N = 12 is around the 
maximum that allows simulations with LRS in rea- 
sonable time (around one hour per instance). For 
each polytope then we have computed the two en- 
tropies Slrs and Sep with both algorithms, fix- 
ing the same maximum value for the discretization 
^max = 1Q24 for all fluxeg ^ 

In Figure |1.2| we show how the quality of the 
BP measure is affected by the discretization, by 
displaying the histogram of the relative differences 
fig _ Sbp-Slrs w jth an increasing number of bins 

&LRS 

per variable g max = 16, 64, 256, 1024. One can see 
how a finer binning of messages improves the qual- 
ity of the approximation, seeminlgy converging to a 



3 



single distribution of errors. It is expected that for 
larger N the histograms would shrink: upon increas- 
ing the number of fluxes, loops become larger and 
the overall topology of the graph becomes more lo- 
cally tree-like, validating the hypothesis behind the 
Bethe approximation. Unfortunately, the huge in- 
crease of computer time experimented in the calcu- 
lation of the volumes using LRS made impossible to 
test systems large enough to make any reasonable 
scaling analysis. 

Finally we address the issue of the computational 
complexity of the algorithm which is a crucial one if 
one is interested in approaching real world metabolic 
networks whose size typically is at least 50 times the 
size of the largest network that can be analyzed with 
exact algorithms. In Figure [O] we display the run- 
ning time of both LRS and BP as a function of the 
number of fluxes N. Interesting, LRS outperforms 
BP up to sizes N ~ 12 where the running time of 
LRS explodes exponentially while BP maintains a 
modest almost-linear behavior. 



Distribution of fluxes in Red Blood Cell 

The algorithm was used to obtain the distribution 
of flux values for each of the reactions of the Red 
Blood Cell metabolism. The maximum allowed val- 
ues for the fluxes, as well as the corresponding sto- 
ichiometric matrix were extracted directly from [8]. 
The network contained 46 reactions and 34 metabo- 
lites. Our distributions appear in Figure |1.3| and 
as can be checked by direct comparison they are 
almost identical to those obtained with the Monte 
Carlo method in with the Figure 5 in [8]. However, 
while the Monte Carlo method appears to be quite 
expensive in computational resources (the authors 
of [8] reported one week of computer computation 
in a Dell Dimension 8200 to obtain their distribu- 
tions) our algorithm converged to the same results 
in a couple of minutes of computation on a similar 
machine. 



Analysis of gene knock-out in E-Coli 

Then, we analyze the distribution function of the av- 
erage fluxes of the metabolic network of the E-Coli. 
The network used contains, in its original format, 
1035 reactions and 626 metabolites, [11]. From this 
network we eliminated all those reactions with frac- 
tional stoichiometric indexes, since our algorithm in 
its present form is unable to deal with them, and all 



those metabolites with connectivity larger than 40. 
These metabolites are cofactors of the metabolism 
and their elimination helps to sparse the matrix im- 
proving the convergence of the algorithm. Finally 
we checked for inconsistencies (i.e. metabolites that 
after the previous transformations are only produced 
or consumed), and we obtain a matrix with 1005 re- 
actions and 560 metabolites. We ran our algorithm 
on this network and computed the average fluxes in 
each reaction. The probability distribution function 



(pdf) of these average values appears in Figure 1.4 



As can be clearly seen the distribution is large and 
may be fitted with a power law with exponent close 
to 1, a result that compares very well with previous 
simulations [12, 13] in real networks and with the 
approach followed by Bianconi and Zecchina [14]. 
Moreover, it most be mentioned that a more careful 
analysis of the data may suggest that the distribu- 
tion of averages fluxes has a richer structure. In 



Figure 1.5 is presented the integrated distribution 
function of the average fluxes of the reactions and a 
clear jump appears for v ss 0.5, and smaller ones for 
v w 0.4 and fi w 0.6. Whether these jumps are just 
due to normal statistical fluctuations (and are cor- 
rectly smeared out in the usual binning process done 
to plot the pdf in log-log scale) or reflect relevant bi- 
ological or structural information of the network is 
not yet known. 

Finally, we concentrated our efforts in the study 
the E-Coli central metabolism [11]. To simplify the 
analysis of the results, all cofactors were removed 
with the exception of ATP and ADP that have been 
thus considered as output and input fluxes respec- 
tively. Under these conditions, the network has 
40 metabolites, 74 internal reactions and 10 exter- 
nal fluxes that were considered irreversible following 
their nominal directions. The maximum fluxes of all 
the reactions were considered equal to 1. 

We then studied the influence of knockouts on 
the volume of the solution space. In each simulation 
the values of the maximum flux was kept constant 
and equal to 1 for all the reactions but the one with 
a knockout. For each knockout we assumed a reduc- 
tion of the maximum permitted flux g max equal to 
the 0.25 of the original value. We assumed that each 
reaction is knockout independently, this, despite the 
fact that it is known that some reactions are associ- 
ated with the same enzyme. 

The distribution of entropy change (AS — Sko~ 
So) of the solution space for the different reaction 
knock-outs turns out to show two possibilities: while 
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most of the reactions have AS > —3 and show little 
and similar impact in the size of the solution space 
a relevant fraction of them, has AS < —3 but with 
a very uniform distribution. To understand if there 
is any connection between these reactions with large 
AS and the structure of the network, we show in 
Figure [F6j the changes in S for the different reaction 
knock-outs. In the x — axis we plot the reactions nu- 
merated and the lines indicate the values of A S for 
each reaction. We annotate also these reactions with 
AS < —3. Most of them are associated with the 
glycolysis pathway, being more specific, with those 
reactions of the glycolysis process that show little 
redundancy in the topology of the network. Other 
reactions, like FUM, ACONT, SUC and SUCCDli 
appear in the Krebb cycle and again show little path- 
way redundancy in the network. 



Finally in Figure 1.7 we display the correlations 
between the changes in the entropy for different re- 
action knock-outs and the most probable values v* 
and the average values v of the fluxes in the wild 
network. As can be seen, two kinds of regimes are 
divided by a clear threshold at v ~ 0.6 : A first one, 
for small fluxes, that show consistent and slowly in- 
creasing correlations between v and AS, and a sec- 
ond one for large fluxes, where the correlations in- 
crease rapidly but with large fluctuations. The pres- 
ence of this threshold can be understood noting that 
reactions belonging to the linear (glycolysis) an cir- 
cular (Krebb cycle) pathways are, in the wild cell, 
fast-fhix reactions, with average values for the fluxes 
larger than 0.5. Therefore once they are knock-out, 
the metabolic capacities of the network are highly 
affected. 



Conclusions 

We proposed a novel algorithm to estimate the 
size and shape of the affine space of a non full- 
dimensional convex poly tope in high dimensions. 
The algorithm was tested in specific benchmark, i.e. 
random diluted stoichiometric matrices at a given 
ratio a = M/N and fixed number of terms differ- 
ent from zero K, in each of the reactions, with re- 
sults that compare very well with those of exact al- 
gorithms. Moreover, we show that while the run- 
ning time of exact algorithms increases more than 
exponentially for already moderate sizes, our algo- 
rithm keeps a polynomial behavior for sizes as large 
as N = 120. The program was run on the Red 



Blood Cell metabolism, showing with less compu- 
tational effort, results that compare very well with 
those previously obtained using Monte Carlo meth- 
ods. Then, we calculate the distribution of the av- 
erage values of the fluxes in the metabolism of the 
E-Coli and present results that are consistent with 
those of the literature. Finally, our program was 
used to study the E-Coli central metabolism, and 
we show that, as expected, reactions with little re- 
dundancy are the ones with more impact in the size 
of the space of the metabolic solutions. Specifically, 
most of the reactions associated with the transfor- 
mation of glucose in pyruvate, belong to this set, as 
well as some reactions in the citric cycle. In addition 
we show strong correlations between the characteris- 
tics of the flux distributions of the wild type network 
and the changes in size of the space of solutions after 
reaction knock-outs. 

Let us conclude by noting that in principle the 
presented approach can be extended to deal with 
constraints whose functional form is more general 
than linear, provided that the number of variables 
involved in each of the constraints remains small, as 
in the case of inequalities enforcing the second low 
of thermodynamics for the considered reactions [15]. 
Work is in progress in this direction. 



Methods 

Volume computation 

From a computational point of view, the problem of 
the exact computation of the volume of a polytope 
with current methods requires the enumeration of 
all its vertexes. The vertex enumeration problem is 
#P-hard [16,17], but even the problem of computing 
the volume, given the set of all vertexes is a big com- 
putational challenge. Various algorithms exist for 
calculating the exact volume of a polytope from its 
vertexes (for a review see [18]), and many software 
packages are available in the Internet. Computa- 
tional limitations restrict however exact algorithmic 
strategies to cope with polytopes in relatively few di- 
mensions (e.g N — M around 10 or so). To overcome 
such severe limitations we will introduce a very effi- 
cient approximate computational strategy that will 
allow us to compute the volume and the shape of the 
space of solutions for real- world metabolic networks. 

Although the coefficient A could be explicitely 
calculated, it turns out that as far as only relative 
volume quantities are concerned, as in the case of the 
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in silico flux knock-outs introduced below, this term 
factors out and therefore we will drop it from the rest 
of the computation. If one is not interested in con- 
stant prefactors, one can proceed as follows for the 
integration of Eq. [6] Consider the regular orthogonal 
grid A e of side e partitioning R . This grid maps via 
into a partition T e of The number of 

cells Af t of A e intersecting II is equal to the numbers 
of cells of r e intersecting </> -1 (II). Finally, the vol- 
ume in Eq.|5]is proportional to lim c _>o e N ~ M Af e - The 
same applies to the computation of the marginal of 
Eq. [7] now noting that the constant pre-factor does 
not depend on v^. The reader may wonder why one 
does not directly integrate Eq. [5] by using its right 
hand side, the problem is that the change of coor- 
dinates induced by <j> may destruct the sparsity of 
the original matrix s a ,i, a condition that turns out 
to be crucial for the kind of approximation that we 
will use. 

When dealing with integer coefficients Si. a , as the 
ones appearing in stoichiometric relations, a further 
simplification in the approximate volume computa- 
tion is possible: one can (always ingoring a constant 
pre-factor) restrict further to integer solutions of the 
system of equations. Geometrically, the space is tiled 
with small hypercubes and we are actually counting 
the number of hypercubes exactly fulfilling the sto- 
ichiometric equations. In summary, for any e the 
computation of an e-approximation of the volume 
has been recast into a discrete combinatorial opti- 
mization problem that can described (with a slight 
abuse of notation) by the same Eqs. [2] with now dis- 
crete variables Vi £ {0, 1, <?™ ax }, for g™ ax equal to 
the integer part of g max x jy™ ax ; where the integer 
g mM is the granularity of the approximation. 

Belief Propagation 

The metabolic problem can be cast into a constraint 
satisfaction framework where the M stoichiometric 
relations impose a constraint onto a subset of the 
metabolic fluxes. Let A be the set of equations and 
I the set of fluxes. Consider the a-th row of S, and 
let {ii,..., i na } = {i £ a} C I be the labels of the 
fluxes involved in the considered equation having 
stoichiometric coefficients different from zero. Let 
also {di , • • • , a n .} = {a £ i} c Abe the labels of the 
equations in which flux i participates. The emerg- 
ing structure is a bipartite graph, with two types 
of nodes: variable nodes representing the fluxes of 
the reactions and factor nodes imposing mass con- 



servation. In this case marginals become g-modal 
probability densities that for large values of qf 1 ax 
will approximate better and better the continuous 
set of probabilities. 

Under the hypothesis that the factor graph is a 
tree it can be shown [19, 20] that a given flux vector 
v satisfying all flux-balance constraints can be ex- 
pressed as a product of flux and reaction marginals 
[20-22]: 

V{v)=T[P a {{n}iea)J[Pi{vi) 1 - di (8) 

a£A iei 

where di is the number of equations in which flux i>i 
participates (i.e. the degree of site i). The marginal 
probabilities are defined as: 

Pi{vi) = £ V(v) 

Pattnhea) = E • ( 9 ) 

Belief Propagation (BP) is a local iterative algo- 
rithm that allows for the computations of marginal 
probability distributions, which is exact on trees, 
and perform reasonably well on locally tree-like 
structures [20-23]. This approximation scheme al- 
lows the computation of the (logarithm of the) num- 
ber of solutions via the entropy that can be expressed 
in terms of flux marginals: 

V 

= E E Pa({Vjea})logP a ({v je a}) 

a<£A{u 3 } jea 

iei i^i 

One may wonder how such an approach could be 
useful in a real-world situation where the graph is 
not a tree. Interestingly enough, metabolic networks 
are sparse, i.e. the number of metabolites that typ- 
ically participate to a certain reaction is small with 
respect to the number of metabolites M, moreover 
one can reasonably assume the typical loop length 
to be large enough to ensure weak statistical depen- 
dence of neighboring sites which lay at the heart of 
the Bethe approximation [24,25]. The algorithm is 
based on two type of messages exchanged from vari- 
able nodes to functional nodes, and vice versa: 

• /ii^ a (v): the probability that flux i takes value 
v in the absence of reaction a. 
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• m a ^i (y) : the non-normalized probability that 
the balance in reaction a is fulfilled given that 
flux i takes value v. 

The two quantities satisfy the following set of func- 
tional equations: 

m a ^i(vi) = ^2 W-»o(^) Sj^sajvi ; b a 

{"ilie.\i ; ea 

PH^aiyi) = C'i^a Y\_ m 6-»«( M i) (11) 
b£i\a 

where Cj_> a is a constant enforcing the normaliza- 
tion of the probability /ij_ >a (i/) and <5(- ; •) is the Kro- 
necker delta function. The set of equations 11 can 



be solved iteratively and upon convergence of the 
algorithm one can compute the marginal flux distri- 
butions as: 

Pa{{vi}iea) = M-afa) 6 (Y1 

{vi}ie a l£a 

H v ) = II TO ^( l/ ) ■ ( 12 ) 

lei 

A brute force integration of the discrete set of equa- 
tion would be much too inefficient for analyzing large 
networks, due to the multiple dimensional sum over 
{0, . . . , 9; max }iea\i m the previous equation. A rel- 
evant speed-up in the convergence of these equa- 
tions can be achieved by noting that the convolution 



product in 11 can be efficiently solved with a recur- 



sion. Let us relabel the set of fluxes in equation a as 
{I e a} = {1, . . . , n a } and let assume for the sake of 
clarity that b a = 0: 



R (k+1} M = E R(k} (" " *«,**)***-«(*) (1 3 ) 



for k £ {Q, ...,n a — 1} and initial condition 
i?(°)(0) = 1. The last step of the recursion gives 
us m a ^i{v) = BS na '(v). The increase of perfor- 
mance is substantial since brute force integration 
involves a sum over Jl;ea\i q ; max terms, while itera- 
tion in Eq. refeq:iteration scales just as J2iea\i <7z max - 
In some cases a faster convergence was met follow- 
ing a refinement strategy in the variable q" lax . In- 
stead of choosing a large gf ax from scratch and to 
solve the BP equations from random initial condi- 
tions, the simulation starts using small g™ ax that 
are increased after convergence. Each time <7f iax is 
increased the BP equations are solved starting from 



a function that fits the previous solution. With this 
strategy, networks as large as 40 metabolites and 
120 reactions could be simulated in a couple min- 
utes using a standard laptop. It should be noted 
that when the number n a is large, the computations 
of m a ^i can be done substantially faster by means 
of discrete Fourier transforms, reducing the compu- 
tation time of all messages m a ^i for i — > a from the 
needed n a x J2iea\i ^ max operations to just around 

2x(E iea 9r x )iog(£ iea 9r x )- 
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1 Figures 

1.1 Figure (13] - Running time 
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Logarithm of the running time vs. N for LRS algorithm, and BP algorithm. Averages were taken over 
5000 realizations for the smaller lattices and 500 N = 12. 
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1.2 Figure [TT2| - Discrepance histogram 




three histograms are for N = 12, M = 4, and K = 3, at different value of g max = 16, 64, 256, 1024. 
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1.3 Figure [TT3| - Distribution of fluxes in Blood Cells 
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Distributions of the flux values for each reaction in the red blood cell network. They are arranged 
following the same sequence as in [8]. 
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1.5 Figure [TT5| - Integrated distribution of average fluxes 
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Integrated distribution of average fluxes (y) of the reactions for the E-Coli metabolism. N = 1005 and 
M = 506. Note the jumps for v = 0.4, v = 0.5 and v = 0.6 
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1.6 Figure [TT6| - Impact of knockout in E-Col 
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Changes of S for different reaction knockouts in the central metabolic network of the E-Coli. Annotated 
reactions belong either to the glycolisis pathway or the Krebb cycle 
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1.7 Figure [T77] - Correlation of impact and value 
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Correlations between the change in entropy AS after reaction knockouts and the average values (P) and 
the most probable values (v*) of the fluxes in the central metabolism of the E-coli. The lines are guides to 
the eyes 
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